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Abstract 

We report a computational study of the adsorption of zoledronic acid molecule on hydroxyapatite (001) surface within 
ab initio density functional theory. The systematic study has been performed, from hydroxyapatite bulk and surface, 
and zoledronic acid molecule to the adsorption of the molecule on the surface. The optimized bond lengths and bond 
angles were obtained and analyzed, giving an evidence of structural similarity between subjects under study. The 
formation energies of hydroxyapatite ( 001 ) surfaces with two kinds of terminations were computed as about 1.2 and 
1.5 J/m^ with detailed atomistic structural information. We determined the adsorption energies of zoledronic acid 
molecule on the surfaces, which are -260 kJ/mol at 0.25 ML and -400 kJ/mol at 0.5 ML. An atomistic insight of 
strong binding affinity of zoledronic acid to the hydroxyapatite surface was given and discussed. 
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1. Introduction 

Bisphosphonates (BPs) are widely used as the most 
powerful drug for protecting bone and treating bone dis¬ 
ease such as osteoporosis and other metabolic bone dis¬ 
orders QllSIill. From the 1960s onward, a great 
number of extensive studies on BPs with respect to the 
development and clinical treatment have been carried 
out, but the research held on BPs is still very active, 
being continued to evolve towards a better understand¬ 
ing of treatment mechanisms and a hnding of more po¬ 
tent BP on the basis of it. In this context, zoledronic 
acid (or zoledronate, ZOD), whose chemical name 
is [l-hydroxy-2-(lH-imidazol-l) ethylidene] bisphos- 
phonic acid or 2 -(imidazol-l-yl)-l-hydroxy-ethane- 1 , 1 - 
diphosphonic acid, was developed as a third-generation 
BP and was approved for the oral treatment of bone dis¬ 
ease in 2012 in DPR Korea after the success of its syn¬ 
thesis in our own way by two of the authors (Yong-Man 
Jang and Song-Un Kim). 

BPs are characterized by two phosphonate groups, 
central carbon atom in between, and two side groups R 1 
and R2, while human bone is composed of complex ar¬ 
ray of hydroxyapatite (HAP, Cag(P 04 ) 30 H) crystallites 
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with nano sizes ranging from 30 to 200 nm embedded 
within the collagen matrix lS0i. It is well estab¬ 
lished that the function of BPs to inhibit bone resorp¬ 
tion is originated with their ability of binding strongly 
to bone mineral, that is, HAP crystal, and stiff resis¬ 
tance to chemical and enzymatic hydrolysis |01. In more 
detailed explanation, the P-C-P backbone of BPs made 
from two phosphonate groups and carbon atom has a 
high binding affinity for the compositions of HAP - 
tetrahedral PO 4 groups, OH groups and Ca ions - with 
an additional contribution from the R1 side group (in 
most cases including ZOD, it is -OH) US. Moreover, 
the P-C-P group of BPs is considerably more resistant to 
the dissolution of HAP crystal than the P-O-P group of 
pyrophosphate, of which BPs are stable structural ana¬ 
logues ESiniElElH. 

Varying the another side group R2 of BPs can result 
in differences in antiresorptive potency with several or¬ 
ders of magnitude. It was observed that more potent BPs 
posses a primary, secondary or tertiary nitrogen atom in 
the R2 side chain IBS. At present, the most potent 
antiresorptive BPs include a heterocyclic R2 side chain 
containing a nitrogen atom like risedronate and zole¬ 
dronate ilTn . The structure-activity relationship anal¬ 
ysis showed that the strong binding affinity of zole¬ 
dronate for HAP is related to its 3D shape and atomic 
orientation, indicating an important role of 3D shape of 
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nitrogen-containing BP and the orientation of its nitro¬ 
gen in binding affinity for HAP |@1 ■ 

Atomistic modeling and simulations are a powerful 
tool to describe the structural characteristics of BPs and 
bone mineral HAP and the interaction between them at 
atomic scale, as proved in materials science and molec¬ 
ular science through a vast number of applications. For 
instance, molecular dynamics (MD) simulations based 
on the well-constructed classical force field have pro¬ 
vided the structural information and energetics of bone 
mineral, HAP crystal and its surfaces 1 18, 2^ 21]. 

Bhowmik et al. MlSll obtained the structural parame¬ 
ters of monoclinic HAP crystal and its surface ener¬ 
getics using the consistent valence force field (CVFF), 
presenting a valuable description of the interaction be¬ 
tween polyacrylic acid and HAP. The surface energet¬ 
ics of HAP crystalline surfaces using ab initio density 
functional theory (DFT) calculations within the gener¬ 
alized gradient approximation (GGA) for the exchange- 
correlation functional has been studied by Zhu and 
Wu d, testing the effects of slab thickness, vacuum 
width between slabs and surface relaxation on surface 
energy. Barrios 1211] has investigated the interaction 
between collagen protein and HAP surface by using 
a combination of computational techniques, DFT and 
classical MD methods. Duarte et al. 12211 performed 
molecular mechanics simulations for molecular struc¬ 
tures of 18 novel hydroxyl- and amino-bisphosphonates 
to examine the interaction between BPs and hydroxya¬ 
patite and to extract relating structural characteristics of 
BPs and their affinities for the mineral, which are in 
agreement with in vitro and in vivo studies for some of 
the studied BPs. To the best our knowledge, however, 
investigation on surface adsorption of zoledronate on 
HAP surface with its detailed atomistic structure based 
on quantum mechanics is still scarce, in spite of such 
extensive theoretical studies of BPs and HAP surface, 
and we believe that ab initio study on this phenomenon 
should definitely contribute to a better understanding 
of the interaction of zoledronate with bone mineral at 
atomic and electronic scale. 

In this paper, we carry out systematic ab initio DFT 
calculations for zoledronate molecule, hydroxyapatite 
bulk and surface, and surface adsorption of zoledronate 
molecule on hydroxyapatite surface. Our special focus 
is placed on the adsorption of zoledronate molecule on 
hydroxyapatite surface, providing the adsorption energy 
and atomistic structures of surface complexes, and an 
insight how charge transferring is occurred in the event 
of adsorption. This is aimed to get a reliable insight 
for bone protection effect of zoledronate at electronic 
scale. In the following, we first describe a computa¬ 


tional method in Sec.|2] present the results for structural 
parameters and electronic properties of hydroxyapatite 
bulk and surface, zoledronate molecule, and for binding 
of zoledronate to hydroxyapatite surface in Sec. [2 and 
finally give our conclusions in Sec.|4] 


2. Computational Method 


For the DFT calculations in this work, we have 
employed SIESTA code ia which solves numerically 
Kohn-Sham equation within DFT |3 25] using a lo¬ 
calized numerical basis set, namely pseudo atomic or¬ 
bitals (PAG), and pseudopotentials for describing the 
interaction between ionic core (nucleus plus core elec¬ 
trons) and valence electrons. The BLYP GGA func¬ 
tional (the Becke exchange functional ll2^ in conjunc¬ 
tion with the Lee-Yang-Parr correlation functional 12711 ') 
was used for exchange-correlation interaction between 
electrons. For all the atoms, Troullier-Martins 128 1 
type norm-conserving pseudopotentials were generated 
within local density approximation (EDA) 1291] . and 
checked carefully. The basis sets used in this work were 
the DZP type (double f plus polarization). The mesh 
size of grid, which is controlled by energy cutoff to set 
the wavelength of the shortest plane wave represented 
on the grid, has taken a value of 200 Ry. Non-fixed 
atoms were allowed to relax until the forces converge 
less than 0.02 e'V/A. 

We first determine the crystal lattice parameters of 
bulk HAP by performing structural relaxation with 
atomic coordinates optimization using the conjugate 
gradient scheme and appropriate Monkhorst-Pack k- 
points. Through a comparison of the lattice constants 
with the experimental values, we have a confidence of 
the above mentioned selection for calculation parame¬ 
ters. Structural optimization for isolated ZOD molecule 
is performed as well, and this for several conformations 
to select the most stable one, of which total energy is 
the lowest among studied conformations. 

We then select the interesting surface index as (001) 
on the basis of experimental evidence that this face 
provides binding sites for many ionic species M 
3 in . In order to simulate the surface within a three- 
dimensional simulation code, we use two-dimensional 
periodic slabs, which have a thickness of atomic layers 
of 2 crystal unit cells plus 30 A vacuum layer, allowing 
the atoms in one half of a crystal unit cell at each side 
of the slab to relax. These structural parameters in the 
supercell can give well-converged result; increasing the 
thickness of slab to 3 and 4 crystal unit cells and the 
vacuum width up to 50 A result in variations of surface 
energy smaller than +0.05 J/m^. To check the stability 
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of the surface, the surface formation energies are calcu¬ 
lated approximately as the difference between the total 
energies of the slab and the corresponding bulk crystal: 

r = [^slab - ^^bulkj /2A, (1) 


where A^siab and A^buik are the numbers of atoms in the 
surface slab and in the bulk unit cell, A is the slab sur¬ 
face area, and iSsiab and fsbuik are the total energies of 
the slab and bulk structures, respectively 0211 . The sign 
of surface formation energy is a test of surface stabil¬ 
ity: positive (negative) means energy should be pro¬ 
vided (released) in order to create a surface. We note 
that there might be several possible cutting planes for a 
particular Miller index while guaranteeing neutrality of 
the surface charge perpendicular to the surface !^ . 

Final step of the work is to simulate the adsorption 
of ZOD molecule on the relevant relaxed HAP (001) 
surfaces. In order to have the rough, initial surface ad¬ 
sorption structure, we use classical molecular mechan¬ 
ics (MM) method. In this MM simulation. General 
Utility Lattice Program (GULP) |34| is utilized. For 
describing the interaction between atoms, we use the 
Dreiding force field, where the potential energy is de¬ 
scribed as the sum of the contributions resulting from 
the bonded interactions (bond stretching, bond bend¬ 
ing, and torsions) and from the non-bonded interactions 
(e.g., electrostatic interaction, and van der Waals inter¬ 
action) (3. After making a guess for the initial state 
of the adsorption complex, we also carry out atomic re¬ 
laxation for the surface complex - ZOD molecule and 
surface atoms in the slab - to obtain a stable final state. 
Then, the adsorption energy can be calculated as fol¬ 
lows: 


^ads “ [^mol-slab (^slab "b A^mol^mol)] 5 ( 2 ) 

A^mol 

where N^oX is the number of adsorbed molecules, and 
iimoi-siab and Zimoi are the total energies of the adsorbed 
surface and of the isolated molecule, respectively. The 
adsorption energy can be either negative or positive: 
negative (positive) means energy should be released 
(provided) during the adsorption, indicating that the ad¬ 
sorption is (not) spontaneous exothermic (endothermic) 
process. 


CajQ(P 04 )g( 0 H) 2 . In fact, there must be four hydroxyl 
groups in the unit cell with the P6^lm space group, 
each o xyg en atom of hydroxyl group with 1/2 occu¬ 
pancy 12111 . To make the DFT calculations enable, there¬ 
fore, we have made a model with full occupancies for 
the hydroxyl groups, by assigning alternate 0 and 1 oc¬ 
cupancies to these hydroxyl groups, which results in the 
change of the space group from P 63 /OT to P 63 . How¬ 
ever, the modified model has a net electric polarization 
contrary to the experiment which shows zero polariza¬ 
tion, because all the hydroxyl groups in the model are 
oriented in the same direction. To mimic the real struc¬ 
ture where exists disorder in the relative orientation of 
the parallel OH channels, so that electric polarizations 
are compensated each other, we have created a super¬ 
cell containing two unit cells in the [ 100 ] direction, and 
assigning opposite orientations to the two OH channels 
in the supercell. The crystal structure with antiparal¬ 
lel hydroxyl groups in a double unit cell compared to 
the original hexagonal structure is monoclinic with P2i 
space group. We have used this structure in this work, 
thus allowing a direct comparison of the simulation re¬ 
sults with experimentally determined surfaces. 




■[ 100 ] 


Figure 1: (Color online) Top view (a) and side view (b) of supercell of 
bulk hydroxyapatite crystal with monoclinic P2i space group, which 
is doubled the unit cell in [100] direction so as to compensate the 
electric polarization by assigning opposite orientations to the two OH 
channels in hexagon surrounded by triangle of Cal ions, (a) top view 
and (b) side view. (Ca: Green, P: purple, O: red, and H: white) 


3. Results and discussion 

3.1. Bulk hydroxyapatite 

The crystal structure of hydroxyapatite is hexagonal 
with space group Pb^Jm, which contains a formula unit 


The full optimization of the structure, allowing the 
cell shape and volume, and ionic positions to relax, was 
performed. We have checked the convergence with re¬ 
spect to the special k-points; increasing the k-points set 
from (1 X 2 X 3) to (2 X 4 X 6 ) leads to the change be- 
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tween the total energies as about 5 meV, thus indicat¬ 
ing the safe use of (1 x 2 x 3) set. Figure [T] depicts the 
fully optimized supercell of bulk HAP in this work. The 
calculated structural parameters and chemical bonding 
properties of bulk HAP are given in Table [1] The lat¬ 
tice constants (a=9.348, b—9352, and c=6.955 A) of the 
fully optimized structure are in good agreement with the 
experimental values {a-b—9A3, c=6.891 A) |36, 37] 


(less than 1.0% error) as well as with the earlier DFT 
results ^ 21). Note that a is the half of lattice 
constant of the supercell. 


Table 1: Structural parameters of bulk hydroxyapatite crystal struc¬ 
ture with monoclinic P2\ space group, compared with experimental 
results._ 


Lattice parameters 


This work 

Experiment" 

a, b, c (A) 
a,/3,jn 

9.348, 9.352, 6.955 
90, 90, 119.86 

9.430, 9.430, 6.891 
90, 90, 120 

Average of bond lengths (A) 

P-0 

1.578 

1.540 

Cal-O 

2.423 

2.478 

Cal-On 

2.339 

2.354 

Ca2-0 

2.553 

2.556 

Average of bond angles (°) 

O-P-0 

109.44 

109.45 

O-Cal-0 

98.62 


0-Ca2-0 

99.40 



° Ref. 


As shown in Figure [T] there are two different Ca sites 
(denoted as Cal and Ca2), four oxygen sites (Ol, 02, 
and 03 of phosphate tetrahedron, and Oh of hydroxyl 
group), one P site, and one H site. The Cal atom is coor¬ 
dinated to seven oxygen ions, which are six oxygens of 
different five PO 4 groups and one oxygen of OH group, 
while the Ca2 atom has a ninefold coordination with 
oxygen ions situated in six different PO 4 tetrahedra. 
The three Cal atoms also form a triangle at the same 
plane normal to the c axis, whose center is occupied 
by an OH group, and the two triangles form a hexag¬ 
onal screw configuration. The calculations yield aver¬ 
age bond lengths as 1.578 A for P-O bond, 2.423 A for 
Cal-O, 2.339 A for Cal-On, and 2.553 A for Ca 2 - 0 , 
compared to the corresponding experimental values of 
1.54, 2.478, 2.354, and 2.556 A, respectively. The sim¬ 
ilar agreement is obtained in the bond angles, as shown 
in Table [1] These comparisons confirm that our com¬ 
putational parameters are reasonably satisfactory for a 
good assessment of surface calculations. 


3.2. Hydroxyapatite (001) Surface 

In this work, we have selected the Miller indices of 
HAP surface as (001), since there exists experimental 
evidence that the ( 001 ) surface is the most stable among 
several surfaces with low indices and provides binding 
sites for many ionic species ii. As we use the super¬ 
cell doubled in the [ 100 ] direction as a unit cell for bulk 
HAP, the (001) surface unit cell is an oblique (2x1) sur¬ 
face cell and there might be two kinds of terminations; 
(1) 2 Ca 2 -( 6 P 04 - 6 CaT 20 H)- 2 Ca 2 , denoted as Type 
(I), and (2) ( 3 P 04 - 3 Cal- 0 H)- 4 Ca 2 -( 3 P 04 - 3 Cal- 0 H), 
denoted as Type (II). They are repeated in the [001] di¬ 
rection, assigned to building layer, and belong to the 
Tasker (Il)-type surface with no electric dipole mo¬ 
ment perpendicular to the surface. 


Table 2: Surface energies (J/m^) according to vacuum thickness and 
building layers of hydroxyapatite (001) surfaces with (2x1) surface 
unit cell. The values in parenthesis represent the number of atoms. 


Type 

(I); 2Ca2- 

( 6 P 04 - 6 CaL 20 H)- 2 Ca 2 


Size 

Unrelaxed 

Relaxed Ref." 


20 

1.460 

1.208 

Vacuum (A) 

30 

1.461 

1.205 


40 

1.448 

1.197 


50 

1.449 

1.205 


4(176) 

1.461 

1.205 

Layer 

6 (264) 

1.474 

1.218 


8 (352) 

1.496 


Type (II); ( 3 P 04 - 3 Cal 0H)-4Ca2- 

■( 3 P 04 - 3 Cal 0 H) 


4(176) 

2.084 

1.514 1.01 

Layer 

6 (264) 

2.057 

1.506 


8 (352) 

2.078 



° Ref. I2II1 


The HAP (001) surface has been modelled using 
three-dimensional periodic supercell (slab) with two 
equivalent surfaces at bottom and top side. To ensure no 
interaction between bottom surface of the above image 
slab and top surface of the present slab across the vac¬ 
uum region, the vacuum region must be wide enough. 
And the atomic layer, which is consisted of surface layer 
allowed to relax and crystal layer fixed at its crystalline 
position, should be thick enough so that the two sur¬ 
faces of each slab do not interact through the crystal 
layer. To check the convergence according to the vac¬ 
uum region, we tested 20 , 30, 40, and 50 A thickness, 
and confirmed that there is no distinct change between 
20 and 50 A thick vacuums (Table | 2 ]i. Therefore, the 
vacuum region of 30 A thickness will be used in the 
following calculations. The thickness of the slab is usu¬ 
ally expressed in terms of a number of building layers. 
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Figure 2: (Color online) Top view (upper panel) and side view (lower panel) of optimized structure of hydroxyapatite (001) surface with 
(2x1) surface unit cell, vacuum thickness of 30 A, 4 building layers, and termination 2Ca2—( 6 P 04 - 6 Cal- 20 H)- 2 Ca 2 (a) and termination 
( 3 P 04 - 3 Cal- 0 H)- 4 Ca 2 —( 3 P 04 - 3 Cal- 0 H) (b). The top and bottom layers are allowed to relax, while the layers between dashed lines are fixed at 
their crystalline positions. (Ca: green, P: purple, O: red, H: white) 


where one layer contains 44 atoms. The four layers (two 
bulk unit cells, 176 atoms), six layers (three unit cells, 
264 atoms), and eight layers (four unit cells, 352 atoms) 
were tested with the vacuum width of 30 A, and the 
change of surface energies between 4 layers slab and 
8 layers slab is only 0.05 J/m^. Thus we will use a slab 
with 4 building layers in the study of surface adsorption. 

As listed in Table ^ the surface energy of Type (II) 
surface is slightly higher than that of Type (I) surface, 
indicating the Type (I) terminated surface is more fa¬ 
vorable than the Type (II) terminated surface. We see 
that the surface relaxation in the Type (I) surface is not 
really much compared to the Type (II) surface, since the 
difference between unrelaxed and relaxed surface ener¬ 
gies in the former case is not remarkable contrary to the 
latter case. Since there is no data available in the litera¬ 
ture for Type (1), we only compared the calculated Type 
(11) surface energy with the previous SIESTA work 12III . 

Figure shows the fully relaxed atomistic structure 
of the HAP (001) surfaces modelled by a slab with vac¬ 
uum thickness of 30 A and four building layers. It is ob¬ 


served that the coordination number (CN) of Cal atom 
changes from 7 in bulk to 6, and CN of Ca2 from 9 to 6 
in the Type (I) termination, while in the Type (II) termi¬ 
nation the CN of Cal atom changes from 7 to 6, and Ca2 
atom from 9 to 5, and 6. The coordination environments 
around the surface oxygen atoms (Ol, 02, 03, and Oh) 
are also changed from the their bulk environment, while 
P atoms are still fully surrounded by four oxygen atoms 
like in bulk. 

Table [3] shows the bond lengths of cation-oxygen and 
bond angles of oxygen-cation-oxygen at the top surface 
layer. In the case of Type (II) surface, we consider the 
two kinds of Cal and three kinds of Ca2 surface atoms. 
In both cases, the averages (1.583, 1.585 A) of bond 
lengths of P-O are a bit expanded compared with the 
bulk (1.578 A), and the averages of bond angles get 
smaller than in the bulk. The bond lengths of Ca2-0 
are clearly contracted at the surface; 2.434 in Type (I), 
and 2.325,2.380, and 2.414 A in Type (II), which are all 
smaller than the bulk value 2.553 A. From these obser¬ 
vations, it can be concluded that the undercoordinated 
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Table 3: Bond lengths of cation-oxygen and bond angles of oxygen- 
cation-oxygen at hydroxyapatite (001) surface with coordination 
numbers (CN) of cations. 


Type (I); 2Ca2-(6PQ4-6Cal-2QH)-2Ca2 



Length (A) 

Angle (°) 


CN 

Range 

Average 

p 

1.583 

103.47-116.23 

109.39 

4 

Cal 

2.352 

64.54-155.26 

101.50 

6 

Ca2 

2.434 

62.45-137.08 

94.01 

6 

Type (II); ( 3 P 04 - 3 Cal-OH)- 4 Ca 2 -( 3 P 04 - 3 Cal-OH) 


P 

1.585 

102.07-117.07 

109.22 

4 

Cal 

2.420 

60.97-159.10 

99.72 

6 

Cal 

2.343 

64.73-160.44 

102.36 

6 

Ca2 

2.325 

66.52-147.60 

101.91 

5 

Ca2 

2.414 

60.91-170.66 

98.61 

6 

Ca2 

2.380 

61.52-154.20 

99.52 

6 


surface Ca atoms and O atoms can be favorable adsorp¬ 
tion sites. 

3.3. Zoledronic acid molecule 

To simulate an isolated molecule, we have used a cu¬ 
bic supercell with lattice constants of a = = c = 50 A, 

which has been proved to be enough to prevent the ar¬ 
tificial interaction between adjacent molecules. There 
might be four different conformations, as presented in 
Figure|3] We can make a distinction between the confor¬ 
mations according to the directions of two pairs of OH 
groups attached to the two P atoms; (1) ZODboui for the 
case where the directions of both OH pairs are outward 
against the nitrogen heterocyclic ring, (2) ZODBin for 
the case of inward directions of both pairs, (3) ZODnih 
for the case of inward direction of one pair in the same 
side of nitrogen atom, and (4) ZODnoui for the case of 
outward direction of one pair in the N side. After the 
total energy minimizations to get the optimized struc¬ 
ture of molecule, we have compared the total energies 
of four conformations. The calculations tell us the most 
stable structure is just the third case, ZODnih, which will 
be used in the following adsorption study. 

The optimized bond lengths and bond angles related 
with P atoms in ZODNin conformation are listed in Ta¬ 
ble |4] The averages of P-O bond lengths are 1.597 
A in PI side and 1.586 A in P2 side, which are simi¬ 
lar to those of bulk HAP (1.578 A) and of HAP (001) 
surface (1.583 A). The averages of O-P-O bond an¬ 
gles (112.11 and 111.42°) are also close to those of bulk 
HAP (109.44°) and HAP surface (109.39°). It is found 
that those values of P2 side are slightly closer to the 
bulk and surface values than PI side, so that P2 side are 




O 



-4885.041 eV 




(1) ZODsout 


0.900 

I_ 


-4885.479 eV 
(2) ZODBin 
0.462 


-4885.941 eV 
(3) ZODNin 


-4885.373 eV 
(4) ZODNout 
0.568 


Figure 3: (Color online) Optimized structures of zoledronic acid with 
four different conformations, (1) ZODboui, (2) ZODgi^, (3) ZODf.ji„, 
and (4) ZODnoui- (C: grey, P: purple, N: blue, 0: red, H: white) 


more favorable to binding with HAP due to the closer 
structural similarity. 

Figure 0] shows the computed isodensity surfaces 
for the highest occupied molecular orbitals (HOMO) 
and lowest unoccupied molecular orbitals (LUMO) of 
ZOD (conformation ZODNin)- While the HOMO and 
HOMO-1 are localized on the heterocyclic ring con¬ 
taining nitrogen atom, LUMO and LUMOh- 1 are on 
the phosphonate groups. This leads to the intramolec¬ 
ular charge separation upon excitation. The change in 
electronic distribution between HOMO and LUMO in¬ 
dicates that two phosphonate groups can play a role of 
electron acceptor in the chemical reaction, while the 
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Table 4: Optimized bond lengths and angles related with phosphorus 
atoms in zoledronic acid conformation, ZOD^jn. Ave. means average. 



P 

related 


Bond length (A) 

Bond angle (°) 

Pl=01 

1.495 

01=Pl-02 

117.26 

Pl-02 

1.654 

01=Pl-03 

115.10 

Pl-03 

1.641 

01=P1-C1 

115.22 

Pl-Cl 

1.902 

02-P1-03 

103.97 



02-Pl-Cl 

102.96 



03-Pl-Cl 

100.12 


Ave. 

O-Pl-O 

112.11 


Ave. 

0-Pl-Cl 

106.10 


P2 

related 


P2=04 

1.495 

04=P2-05 

120.01 

P2-05 

1.622 

04=P2-06 

117.35 

P2-06 

1.641 

04=P2-C1 

99.93 

P2-C1 

1.934 

05-P2-06 

96.91 



05-P2-C1 

112.73 



06-P2-C1 

110.40 


Ave. 

0-P2-0 

111.42 


Ave. 

0-P2-C1 

107.69 


P1-C1-P2 

98.65 


Figure 4: (Color online) Moleculai* orbitals of zoledronic acid confor¬ 
mation, ZODjsji,^. 


heterocyclic ring containing nitrogen atom could be a 
donor. 


3.4. Adsorption of zoledronic acid on hydroxyapatite 
(001) surface 

To begin with adsorption, we have utilized GULP to 
obtain rough structure of ZOD absorbed on HAP (001) 
surface, where Dreiding forcefield was adopted. Simu¬ 
lated annealing was performed, increasing the tempera¬ 
ture from 300 K to 10000 K and subsequently decreas¬ 
ing back to 300 K with the temperature interval of 50 K. 
Using the obtained rough structure as the starting struc¬ 
ture, we then performed atomic relaxation with SIESTA 
code to get the final optimized structure of zoledronic 
acid adsorbed HAP (001) surface. 

In this work we have tried to simulate the adsorption 
of ZOD molecule at 0.5 ML (one molecule on (2x1) 
surface cell) and at 0.25 ML coverages (one molecule 
on (2x2) surface cell). As mentioned above, there are 
two possible terminations in the HAP (001) surface, and 
therefore four kinds of simulation tasks were carried 
out. We have used the slab model with 4 building layers 
and vacuum thickness of 30 A, which guarantee to give 
a reliable adsorption energy. To make a systematic er¬ 
ror small, the pristine surface energies for (2x2) surface 
slabs were also calculated. In Table |5] the adsorption 

Table 5: Calculated adsorption energies (kJ/mol) of zoledronic acid 
on hydroxyapatite (001) surface. 



Coverage 


0.5 ML (2x1) 

0.25 ML (2x2) 

Type (I) 

-388.44 

-246.70 

Type (II) 

-439.49 

-268.34 


energies are listed. All the adsorption energies are neg¬ 
ative, indicating that all the adsorptions are exothermic 
reactions. We see that the adsorption on Type (II) sur¬ 
face is slightly more favorable than on Type (I) at both 
0.5 ML and 0.25 ML, although the Type (I) surface for¬ 
mation from the bulk needs more energy than the Type 
(II) surface. 

Figure |5] shows the optimized atomistic structure of 
ZOD adsorption complexes on HAP (001) surfaces at 
0.25 ML coverage. It is observed that in the Type 
(I) surface the hydrogen atom of ZOD’s phosphonate 
group moved to oxygen atom of HAP surface’s phos¬ 
phate group, making formation of additional OH group 
on the surface, and thus indicating that the adsorption 
is kind of chemisorption caused by proton exchange. 
There are several hydrogen bonds between ZOD’s phos¬ 
phonate OH group and oxygen atoms of the surface’s 
phosphate groups, and vice versa. It is also found that 
Cal and Ca2 atoms with 6 coordination number (CN) 
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Figure 5: (Color online) Top view (upper panel) and side view (lower panel) of optimized atomistic structure of zoledronic acid adsorption 
complexes on hydroxyapatite (001) surfaces with Type (I) (a) and Type (II) (b) terminations at 0.25 ML coverage. (Ca: green, P: purple, O: red, N: 
blue, H: white) 


make bonds with oxygen and nitrogen atoms of ZOD, 
respectively. Such bonding was expected based on the 
surface relaxation analysis and HOMO-LUMO distri¬ 
bution of isolated ZOD molecule. In the case of Type 
(II) surface, there are also hydrogen bonds between the 
surface and ZOD, but the move of hydrogen atom is not 
observed, which indicates that the adsorption is kind of 
physisorption. Nevertheless, the magnitude of adsorp¬ 
tion energy on Type (II) surface is larger than that on 
Type (I) surface. 

To make it clear the charge transfer occurred during 
the adsorption, we show the electron density difference, 
Psuif+moi - (Psurf+Pmoi), plots in Figure^ where (a) is for 
the isosurface figure with the value of +0.04 eV/A^, and 
(b) and (c) are the contour plots on the planes around hy¬ 
drogen bond and Ca cation, respectively. These figures 
clearly give the evidence of charge transferring at the 
event of adsorption, from hydrogen of phosphate group 
of ZOD molecule to oxygen of PH 4 group of HAP sur¬ 
face making formation of hydrogen bond, and from Ca 
of HAP surface to oxygens of phosphonate group form¬ 
ing the coordinate bond. Therefore the adsorption of 
ZOD on the HAP (100) surface is surely chemisorption. 


We calculated the density of states (DOS) of elec¬ 
trons, total and partial DOS shown in Figure|2]and atom 
projected and partial DOS in Figure 0 In Figure |7] we 
see the energy spectrum of isolated ZOD molecule (a), 
the DOS of ZOD molecule that was adsorbed on HAP 
surface (b), the DOS of HAP (100) surface (c), and the 
whole complex system comprised of the adsorbed ZOD 
molecule and HAP (100) surface. Some hybridization 
between ZOD and HAP surface electrons is shown in 
the figures. Which atoms cause the charge transferring 
at the adsorption? To make an answer to this question, 
we carefully investigate the atom projected partial DOS. 
The Is peak of hydrogen of isolated ZOD molecule 
placed over 0 eV (red line in Figure [T] (a) panel) is re¬ 
markably weakened after the adorption (blue line in (a) 
panel), indicating the loss of electrons from hydrogen 
and the gain by oxygen, and thus the formation of hy¬ 
drogen bond between hydrogen atom of ZOD and oxy¬ 
gen of HAP surface. We also see the hybridizations be¬ 
tween hydrogen Is state and oxygen 2 p state in (a) and 
(b) panels, and between oxygen 2s state and calcium 3p 
state in (c) and (d) panels. 
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Figure 6: (Color online) Electronic charge density difference of zoledronic acid adsorption complex on hydroxyapatite (001) surface, (a) isosurface 
figure with the value of ±0.04 eV/A^ (yellow for + and blue for - value), (b) contour figure on the plane containing hydrogen bond between O of 
phosphonate group of HAP surface and hydrogen of phosphate group, and (c) contour figure around Ca cation. (Ca: green, P: purple, O: red, N: 
blue, H: white) 


4. Conclusion 

In conclusion, we have attempted to make a model¬ 
ing of adsorption of zoledronic acid on hydroxyapatite 
(001) surface to get an atomistic insight of bone protec¬ 
tion. The systematic study has been performed, from 
hydroxyapatite bulk and surface, and zoledronic acid to 
adsorption of the molecule on the (001) surface. We 
have carried out the structural optimizations and atomic 
relaxations of the bulk, molecule, surface, and adsorp¬ 
tion complexes on the surfaces, and obtained the struc¬ 
tural informations and energetics. 

Using the three-dimensional periodic supercell 
model, we determined the stable conformation of the 
molecule and calculated the molecular orbitals. It was 
concluded that two phosphonate groups can play a role 
of electron acceptor in the chemical reaction, while the 
heterocyclic ring containing nitrogen atom can be a 
electron donor. After verifying the validity of compu¬ 
tational parameters of SIESTA work through the bulk 


hydroxyapatite calculation, surface modeling and relax¬ 
ations were performed, and surface formation energies 
were calculated for two kinds of (001) surface termina¬ 
tions, which are about 1.2 and 1.5 J/m^. Subsequently, 
the adsorption of zoledronic acid on the relaxed sur¬ 
face was studied, obtaining the surface binding structure 
and calculating the adsorption energies. We found that 
the surface Ca atoms and oxygen atoms of phosphate 
groups can form surface bond including hydrogen bond 
and coordinate bond with nitrogen, hydrogen, and oxy¬ 
gen atoms of zoledronic acid. The calculated adsorp¬ 
tion energies are about -260 kJ/mol at 0.25 ML cover¬ 
age and -400 kJ/mol at 0.5 ML coverage, indicating the 
strong binding affinity of zoledronic acid to hydroxyap¬ 
atite surface. 

We have made an interpretation of such strong bind¬ 
ing affinity through the analysis of atomistic structures, 
electron density difference, the density of states and 
Hirshfeld charges of atoms relevant to surface binding. 
The results showed that bond lengths and bond angles 
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Figure 7: (Color online) Total and partial density of states in iso¬ 
lated zoledronic acid (a), adsorbed zoledronic acid (b), hydroxyap¬ 
atite (100) surface (c) and zoledronic acid adsorption complex on the 
surface (d). 


of phosphonate group of zoledronic acid are similar to 
those of phosphate group of hydroxyapatite bulk and 
surface, indicating the strong binding affinity related 
to the structural similarity. It was also found that the 
charge transferring is occurred mainly on the side of 
phosphonate group in one kind surface, while in another 
surface it is occurred on both imidazol ring contain¬ 
ing nitrogen atom and phosphonate group of zoledronic 
acid. 
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Figure 8: (Color online) Atom projected partial density of states; (a) 
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surface. 


Notes 

The authors declare that they have no conflict of in¬ 
terest. 


References 

[1] F. H. Ebetino, A.-M. L. Hogan, S. Sun, M. K. Tsoumpra, 
X. Duan, J. T. Triftt, A. A. Kwaasi, J. E. Dunford, B. L. Baimett, 
U. Oppermann, M. W. Lundy, A. Boyde, B. A. Kashemirov, 
C. E. McKenna, R. G. G. Russell, The relationship between the 
chemistry and biological activity of the bisphosphonates. Bone 
49 (2011) 20-33. 

[2] R. Bartl, B. Frisch, E. von Tresckow, C. Bartl, Bisphosphonates 
in Medical Practice, Springer-Verlag, Berlin, Germany, 2007. 

[3] F. P. Coxon, K. Thompson, M. J. Rogers, Recent advances in un¬ 
derstanding the mechanism of action of bisphosphonates, Curr. 
Opin. Pharmacol. 6 (2006) 307-312. 


rea. 


10 

































[4] R. G. G. Russell, N. B. Watts, F. H. Ebetino, M. J. Rogers, 
Mechanisms of action of bisphosphonates: similarities and dif¬ 
ferences and their potential influence on clinical efficacy, Osteo- 
poros Int. 19 (2008) 733-759. 

[5] R. G. G. Russell, Bisphosphonates: The first 40 years. Bone 49 
(2011) 2-19. 

[6] T. Hassenkam, G. E. Fantner, J. A. Cutroni, J. C. Weaver, D. E. 
Morse, R K. Hansma, High-resolution AFM imaging of intact 
and fractured trabeculai' bone. Bone 35 (2004) 4. 

[7] P. Alberius-Henning, E. Adolfsson, J. Grins, Triclinic oxy- 
hydroxyapatite, J. Mater. Sci. 36 (2001) 663. 

[8] J. C. Elliot, Structure and Chemistry of the Apatites and Other 
Calcium Orthophosphates, Elsevier Science, Amsterdam, 1994. 

[9] R. G. G. Russell, R. C. Muhlbauer, S. Bisaz, D. A. Williams, 
H. Fleisch, The influence of pyrophosphate, condensed phos¬ 
phates, phosphonates and other phosphate compounds on the 
dissolution of hydroxyapatite in vitro and on bone resorption 
induced by parathyroid hormone in tissue culture and in thy- 
ropai'athyroidectomised rats, Calcif Tissue Res. 6 (1970) 183. 

[10] E. G. Kontecka, J. Jezierska, M. Lecouvey, Y. Leroux, H. Ko- 
zlowski, Bisphosphonate chelating agents. Coordination abil¬ 
ity of 1-phenyl-1-hydroxymethylene bisphosphonate towards 
Cu(2+) ions, J. Inorg. Biochem. 89 (2002) 13. 

[11] H. Fleich, Bisphosphonates: mechanisms of action, Endocr. 
Rev. 19(1998) 80. 

[12] M. J. Rogers, S. Gordon, H. L. Benford, F. P. Coxon, S. P. Luck- 
man, J. Monkkonen, J. C. Frith, Cellular and molecular mech¬ 
anisms of action of bisphosphonates. Cancer Suppl. 8 (2000) 
2961. 

[13] A. J. Roelofs, K. Thompson, S. Gordon, M. J. Rogers, Molec¬ 
ular mechanisms of action of bisphosphonates: current status, 
Clin. Cancer Res. 12 (2006) 6222. 

[14] J. E. Dunford, K. Thompson, F. P. Coxon, S. P. Luckman, F. M. 
Hahn, C. D. Poulter, F. H. Ebetino, M. J. Rogers, Structure- 
Activity Relationships for Inhibition of Famesyl Diphosphate 
Synthase in Vitro and Inhibition of Bone Resorption in Vivo by 
Nitrogen-Containing Bisphosphonates, J. Pharm. Exp. Therap. 
296 (2001) 235-242. 

[15] G. H. Nancollas, R. Tang, R. J. Phipps, Z. Henneman, S. Guide, 
W. Wu, A. Mangood, R. G. G. Russell, F. H. Ebetino, Novel 
insights into actions of bisphosphonates on bone: Differences in 
interactions with hydroxyapatite. Bone 38 (2006) 617. 

[16] M. A. Lawson, Z. Xia, B. L. Barnett, J. T. Triffitt, R. J. Phipps, 
J. E. Dunford, R. M. Locklin, F. H. Ebetino, R. G. G. Russell, 
Differences Between Bisphosphonates in Binding Affinities for 
Hydroxyapatite, J. Biomed. Mater. Res. B 92B (2010) 149. 

[17] J. R. Green, K. Muller, K. A. Jaeggi, Preclinical pharmacology 
of CGP 42 446, a new, potent, heterocyclic bisphosphonate com¬ 
pound, J. Bone Miner. Res. 9 (1994) 745-751. 

[18] R. Bhowmik, K. S. Katti, D. Katti., Molecular dynamics simu¬ 
lation of hydroxyapatite-polyacrylic acid interfaces. Polymer 48 
(2007) 664. 

[19] W. Zhu, P. Wu, Surface energetics of hydroxyapatite: a DFT 
study, Chem. Phys. Lett. 396 (2004) 38. 

[20] K. Matsunaga, A. Kuwabai'a, First-principles study of vacancy 
formation in hydroxyapatite, Phys. Rev. B 75 (2007) 014102. 

[21] N. A. Barrios, A Computational Investigation of the Interaction 
of the Collagen Molecule with Hydroxyapatite, Ph.D. thesis. 
Department of Chemistry, University College London, 2010. 

[22] L. F. Duarte, F. C. Teixeira, R. Fausto, Molecular modeling of 
the interaction of novel hydroxyl- and aminobisphosphonates 
with hydroxyapatite, ARKIVOC (v) (2010) 117. 

[23] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, 
P. Ordejon, D. Sanchez-Portal, The Siesta method for ab ini¬ 
tio order-N materials simulation, J. Phys:Condens. Matter 14 


(2002) 2745. 

[24] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Phys. 
Rev. 136 (1964) B86^B871. 

[25] W. Kohn, L. J. Sham, Self-consistent equations including ex¬ 
change and correlation effects, Phys. Rev. 140 (1965) A1133- 
A1138. 

[26] A. D. Becke, Density-functional exchange-energy approxima¬ 
tion with correct asymptotic behavior, Phys. Rev. A 38 (1988) 
3098-3100. 

[27] C. Lee, W. Yang, R. G. Parr, Development of the Colle-Salvetti 
correlation-energy formula into a functional of the electron den¬ 
sity, Phys. Rev. B 37 (1988) 786. 

[28] N. Troullier, J. L. Martins, Efficient pseudopotentials for plane- 
wave calculations, Phys. Rev. B 43 (1991) 1993-2006. 

[29] J. P. Perdew, A. Zunger, Self-interaction coiTection to density- 
functional approximations for many-electron systems, Phys. 
Rev. B 23 (1981) 5048. 

[30] N. Almora-Barrios, K. Austen, N. de Leeuw, A Density Func¬ 
tional Theory study of the binding of glycine, proline and hy- 
droxyproline to the hydroxyapatite (0001) and (0110) surfaces, 
Langmuir 25 (2009) 5018-5025. 

[31] A. Wierzbicki, H. S. Cheung, Molecular modeling of inhibition 
of hydroxyapatite by phosphocitrate, J. Mol. Stinc-Theochem. 
529 (2000) 73-82. 

[32] C.-J. Yu, J. Kundin, S. Cottenier, H. Emmerich, Ab initio mod¬ 
eling of glass corrosion: hydroxylation and chemisorption of 
oxalic acid at diopside and akermanite surfaces, Acta Mater. 57 
(2009) 5303-5313. 

[33] P. W. Tasker, The stability of ionic crystal surfaces, J. Phys. C 
12(1979) 4977^984. 

[34] J. D. Gale, A. L. Rohl, The general utility lattice program 
(GULP), Mol. Sim. 29 (2003) 291-341. 

[35] A. Leach, Moleculai* Modelling: Principles and Applications, 
Prentice Hall, 2001. 

[36] J. Y. Kim, R. R. Fenton, B. A. Hunter, B. J. Kennedy, Powder 
diffraction studies of synthetic calcium and lead apatites, Aust. 
J. Chem. 53 (2000) 679. 

[37] A. S. Posner, A. Perloff, A. F. Diorio, Refinement of the hydrox¬ 
yapatite structure, ActaCryst. 11 (1958) 308. 

[38] N. H. de Leeuw, Density functional theory calculations of local 
ordering of hydroxy groups and fluoride ions in hydroxyapatite, 
Phys. Chem. Chem. Phys. 4 (2002) 3865-3871. 

[39] D. E. Ellis, J. Terra, O. Warschkow, M. Jiang, G. B. Gonzalez, 
J. S. Okasinski, M. J. Bedzyk, A. M. Rossi, J. G. Eon, A the¬ 
oretical and experimental study of lead substitution in calcium 
hydroxyapatite, Phys. Chem. Chem. Phys. 8 (2006) 967-976. 

[40] N. H. de Leeuw, A computer modelling study of the uptake 
and segregation of fluoride ions at the hydrated hydroxyap¬ 
atite (0001) surface: introducing a Caio(P04)6(OH)2 potential 
model, Phys. Chem. Chem. Phys. 6 (2004) 1860-1866. 


11 


